
gen sample = 1 if e(sample) == 1 // Drop all unused observations. This will speed up the margins command significantly
keep if sample == 1
sum advice if sample == 1 // The following generates data for underlying histogram showing # of obs
local min = `r(min)'
local max = `r(max)'
twoway__histogram_gen advice if sample == 1, freq gen(h x) start(`min') width(1)
replace x = . if h == .
sum advice if sample == 1
global obs = `r(N)'
foreach var of varlist h {
	replace `var' = `var' / $obs * 100
	}

estimates restore $model // The following estimates marginal effects at various levels of X
margins, dydx(treat) at(advice = (`min'(10)`max'))
matrix margins = r(table)' // Saves marginal effects as variables
svmat margins, names(col)
matrix at = r(at)
svmat at, names(at)
gen n = _n
sum n if at3 != .
local max `r(max)'
forvalues i = 1/`max' {
	sum at3 if n == `i'
	replace at3 = `r(mean)' if n == `i' + `max' & at3 == .
	} 

gen zero = 0 // The following code is for the graph
sum ll
local min = `r(min)'
sum ul 
local max = `r(max)'
if abs(`min') > `max' {
	local yrange1 = abs(`min')
		}
		else {
			local yrange1 = `max'
		}
sum h
local yrange2 = `r(max)' * 4


* Graph
graph twoway bar h x if x < 0, barw(1) lc(red) bc(red) fintensity(100) yaxis(2) ///
	||	bar h x if inrange(x,0,40), barw(1) lc(gold) fintensity(100) bc(gold) yaxis(2) ///
	||	bar h x if x > 40, barw(1) lc(green) bc(green) fintensity(100) yaxis(2) ///
    ||  rarea ul ll at3 if se !=., color(gs14) fintensity(100) ///
	||  line b at3 if se != ., clcolor(black) clpattern(solid) lwidth(medthick) ///
	||  line zero at3 if se != ., clcolor(black) clpattern(solid) lwidth(thin) ///
		, xlabel($xlabel, nogrid labsize(3.5)) ///
		ylabel($ylabel, axis(1) nogrid labsize(3.5) angle(horizontal)) ///
        ylabel(, axis(2) nogrid noticks nolabels) ///
        yscale(range(-`yrange1' `yrange1') noline alt) yscale(range(0 `yrange2') noline alt axis(2)) xscale(noline) legend(off) ///
        yline(0, lcolor(black) lwidth(thin)) ///
        xtitle("VAA congruence score", margin(top) size(4)) ///
        ytitle("Treatment effect", size(4)) ///
        ytitle("" , axis(2) size(4) margin(zero)) ///
        ysca(axis(2) titlegap(4)) ///
		xsize(6) ysize(4) ///
		title("$title", margin(bottom)) ///
        scale(*1.3) scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white))

